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We discuss two methods of an exact stochastic representation of the non-Markovian quantum 
dynamics of open systems. The first method employs a pair of stochastic product vectors in the 
total system's state space, while the second method uses a pair of state vectors in the open system's 
state space and a random operator acting on the state space of the environment. Both techniques 
lead to an exact solution of the von Neumann equation for the density matrix of the total system. 
Employing a spin star model describing a central spin coupled to bath of surrounding spins, we 
perform Monte Carlo simulations for both variants of the stochastic dynamics. In addition, we 
derive analytical expression for the expectation values of the stochastic dynamics to obtain the 
exact solution for the density matrix of the central spin. 
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I. INTRODUCTION 

The Markovian dynamics of an open quantum system 
S which is coupled to an environment £^ |l| is convention- 
ally described by a master equation for the open system's 
density matrix ps(t) with a generator in Lindblad form 
It is a well-known feature Hi, HSU of this type 
of master equations that it yields a stochastic represen- 
tation for ps(t) in the form of an expectation value over 
an ensemble of pure state vectors. This means that ps(t) 
can be expressed by 



ps{t)^E{\m)m)\), 



(1) 



where \4'{t)) is a stochastic state vector in the Hilbert 
space of the open system and E denotes the expectation 
value. The great advantage of the stochastic representa- 
tion consists in the fact that it leads to efficient Monte 
Carlo techniques in which one propagates an ensemble of 
pure state vectors in the open system's Hilbert space and 
estimates the reduced density matrix through an appro- 
priate ensemble average. 

Recently, an exact stochastic treatment of non- 
Markovian quantum dynamics has been proposed [loj . 
This method is based on a representation of the density 
matrix p{t) of the total system through an expectation 
value of the form 



p(i) = E(|cI>i(t))($2(i)|). 



(2) 



By contrast to the conventional approach one uses in 
this method a pair of random product vectors = 
tpi{t) 'S) Xi{t) and |$2(t)) ^ tp2{t) (S) X2{t) of the total sys- 
tem. These product vectors follow independent stochas- 
tic time-evolution equations that can be constructed in 
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such a way that the average over the probabilistic dynam- 
ics reproduces the exact Schrodinger or von Neumann dy- 
namics of the total system. As demonstrated in j9^,ilO|| the 
evolution equations for the product vectors \^i{t)) and 
|$2(^)) can be chosen to be relatively simple time- local 
stochastic differential equations, describing a piecewise 
deterministic process or a diffusion process (Brownian 
motion) in Hilbert space. This means that it is possi- 
ble to design a representation of non-Markovian quan- 
tum dynamics involving strong memory effects through 
a Markovian unravelling by means of a pair of indepen- 
dent stochastic product vectors. 

This method bears several advantages. The stochas- 
tic differential equations for the product vectors describe 
a Markovian random process for which efficient numeri- 
cal simulation algorithms are known (see, e. g., Ref. [l| 
and references therein). Since the method is based on 
a direct stochastic representation of the full Schrodinger 
dynamics of the total system it does not rely on the con- 
struction of an approximate effective master equation for 
the reduced density matrix; it does not even require the 
existence of such an equation. Furthermore, the method 
allows, at least in principle, the treatment of arbitrary 
correlations in the initial state. This follows from the fact 
that any initial state p(0) can be represented in the form 
p{0) = E(|$i(0))($2(0)|), where |$i(0)) and |$2(0)) are 
random product states |9|. In particular, the method 
does not presuppose that system and environment are 
initially in an uncorrelated tensor product state. Finally, 
the technique not only allows the determination of the 
reduced density matrix ps (t) but also of multitime quan- 
tum correlation functions of the open system. 

There is an important limitation in the applicability 
of the Monte Carlo algorithms based on the stochastic 
representation ^ which is due to the behavior of the 
statistical fluctuations. As shown in Ref. 9] the fluctua- 
tions of the process and, hence, also the statistical errors 
of the Monte Carlo simulation may eventually grow ex- 
ponentially with time. Therefore, the method can gener- 
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ally be expected to be feasible only for short time scales. 
However, it must be noted that the stochastic dynamics 
of the product vectors and |$2(i)) is by no means 

unique, i. e., there exists an infinite number of stochastic 
evolution equations for which the expectation value ([2]) 
exactly represents the full system dynamics. Recently, 
this freedom in the choice of an appropriate stochastic 
dynamics has been employed to develop optimized Monte 
Carlo algorithms which lead to a drastic reduction of the 
size of statistical errors [lH . 

The structure of Eq. ([2]) is not the only possibility of 
obtaining an exact stochastic representation for the total 
density matrix. In fact, one can construct many other 
random functionals whose expectation values lead to the 
desired equation of motion. Here, we examine an al- 
ternative stochastic formulation which employs a pair of 
random state vectors ^-^^ \fp2{t)) in the open sys- 

tem's Hilbert space and a random operator Rsit) on the 
state space of the environment: 

p{t) ^n\Mt))mt)\ REit)). (3) 

Again, one can construct an appropriate stochastic dy- 
namics for the state vectors and \ip2{t)) and for 
the environmental operator Rsit) that guarantees that 
the expectation value ([3]) exactly satisfies the von Neu- 
mann equation of the total system. 

We start our considerations in Sec. |TT] with a descrip- 
tion of the general concepts underlying the stochastic rep- 
resentations given by Eqs. ^ and ([3]). The correspond- 
ing Monte Carlo simulation techniques will be illustrated 
in Sec. IIIII with the help of a spin star model, i. e., a 
model of a central spin that is coupled to a bath of sur- 
rounding spins [12j . In addition to performing numerical 
simulations, we derive analytical expressions for the ex- 
pectation values ^ and ^ and relate these directly to 
the solution of the Schrodinger equation for the total sys- 
tem. It turns out that the method described by Eq. ([3]) 
is particularly useful for the simulation of the dynamics 
for an infinite number of bath spins. Some conclusions 
are drawn in Sec. IIVI 

II. STOCHASTIC REPRESENTATIONS OF 
NON-MARKOVIAN QUANTUM DYNAMICS 

A. General theory 

We consider an open quantum system with Hilbert 
space Tis coupled to an environment with Hilbert space 
Ti-E- The state space of the composite quantum system 
is given by the tensor product space Hs He- Em- 
ploying the interaction picture we write the Hamiltonian 
describing the system-environment coupling as follows, 

(4) 

a 

The Aa{t) and the Ba{t) are interaction picture opera- 
tors acting in Tis and Ti.E^ respectively. The correspond- 



ing von Neumann equation for the density matrix p{t) of 
the composite quantum system is given by 

J^pit)=-^mt),pit)], (5) 

where we set h = 1. 

Our aim is to construct a stochastic representation of 
the total density matrix p{t) in terms of the expectation 
value 

p{t) - E{R{t)). (6) 

Here, R(t) represents a random operator on the state 
space Tis (^He oi the total system. The stochastic pro- 
cess governing the dynamics of this operator must be 
constructed in such a way that the expectation value ([6]) 
satisfies the von Neumann equation ([S]). It turns out that 
there are many possibilities of constructing a stochastic 
representation which meets this requirement. Of course, 
we do not seek just any stochastic formulation, but the 
intention is to find a stochastic process which leads to 
a considerable simplification of the representation of the 
reduced density matrix 

Ps{t)^ti-Ep{t)=EityER{t)) (7) 

of the open system and which allows an efficient numeri- 
cal implementation of its time evolution (tr e denotes the 
partial trace over He)- In the following we discuss two 
(of many other) such possibilities, in which R(t) follows 
a piecewise deterministic process (PDP) 

B. Stochastic process of the form R — |$i)(<I>2| 

The first possibility of a stochastic representation in 
terms of a PDP is given by taking R{t) to be of the form 

i?(t) = |ci>iW)($2(t)|, (8) 

such that we have: 

p(i) = E(|<i>i(0)(<i>2(i)|). (9) 

|$i(t)) and I $2(0) represent a pair of stochastic state 
vectors of the composite quantum system which are cho- 
sen as direct products of system states ipv{t) and envi- 
ronmental states Xi'(^)' 

\<^u{t)) ^Mt)®Xu{t). u=l,2. (10) 

In view of Eqs. ([9]) and (flO)) the reduced density matrix 
ps{t) [see Eq. d?])] can be expressed in terms of the ex- 
pectation value: 

P5(i)=E(|V.l(i))(^2(OI(X2(i)|Xl(t))). (11) 

By contrast to the standard stochastic unravelling of the 
dynamics of open quantum systems, this representation 
employs an average over the product of two quantities. 
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namely the dyadic \tpi){ip2\ formed by a pair tpi, tfj2 of 
state vectors of the open system, and the scalar product 
(X2|xi) of a corresponding pair of environment states. 

According to the ansatz pUj) the states of the 

total system are direct products at any time t, which 
greatly simplifies the representation of the states of the 
system and the simulation of its dynamics. Of course, the 
exact states are generally entangled. This shows that 
the dynamics of the cannot be described by a 

deterministic time evolution. However, as demonstrated 
in [l^ it is possible to reproduce the dynamics of the 
total density matrix with the help of a random Markov 
process. An appropriate system of stochastic differential 
equations for the ij^,^ and the Xv is given by 



\\Ao.im.\\' 



^Ac^it)- l]i,,dN^,{t) (12) 



and 



E 



\Bc.{t)xA 



(13) 



where / denotes the unit operator. The quantities dNa 
are random Poisson increments which satisfy 



^{dN^,{t))^T^,{t)dt 



and 



dN^^{t)dNi3^{t) = 6^p5,,f,dN^^{t). 
The corresponding rates are given by 

\\A^{t)i^,\\-\\B^{t)xA\ 

11^.11 -llx.!! ' 

and we have defined the total rates 

r,(t) = ^r„,(t). 



(14) 
(15) 

(16) 

(17) 



In view of Eq. the stochastic increments dN^i, (t) 
take on the possible values or 1. According to Eq. (fT4|) 
the case dNav{t) = 1 occurs with probability Tavdt. Un- 
der the condition that dNav{t) = 1 for a particular a and 
V, the other increments vanish, and Eqs. p2)) and ([M]) 
imply that for this particular a and v the state vectors 
Tpi, and Xv perform the instantaneous jumps 



\\Aa-lpt,\\ 



llx.ll 

\'BaXv\ 



■B^Xu. (18) 



Note that these jumps preserve the norm of the state 
vectors. Under the condition that all Poisson increments 
vanish, that is dNai/{t) = for all a and u, we have 
dipi^it) — and dxv{t) — ^vXi^dt. This means that 
remains unchanged during dt^ while Xv follows a linear 
drift. 



Summarizing, ipiyit) is a pure, norm-conserving jump 
process, whereas X;y(i) is a PDP with norm-conserving 
jumps and a linear drift. It is demonstrated in Q that 
any initial density matrix of the total system can be rep- 
resented in the form p(0) ^ E(|$i(0))($2(0)|), and that 
the expectation value ([9]) exactly satisfies the von Neu- 
mann equation ([5]). These facts enable us to simulate the 
full non-Markovian quantum dynamics through a Monte 
Carlo simulation of the stochastic differential equations 
Uni) and (HID. 



C. Stochastic process of the form 7? = \iJi){tp2\ ® Re 

The second possibility of a stochastic representation is 
obtained if we take R{t) to be of the form 

R{t) = \,Pi{t)){Mt)\®RE{t), (19) 

such that we have: 

Pit) ^T^i\Mt)){Mt)\(^ REit)). (20) 

In this case we represent R{t) through a pair ipi, -02 of 
state vectors of the open system and a random operator 
Re on the state space of the environment. The reduced 
density matrix of the open system can thus be written as 

psit) ^E{\ij,{t)){Mm^ERE{t)) . (21) 

An appropriate system of stochastic differential equation 
which reproduces the exact von Neumann dynamics for 
the expectation value (f20|) is given by 

dip,{t) = ^{-iLo^^A^- I)i^,{t)dN^,{t), (22) 

a 

dRE{t) = VRE{t)dt 

+ (MaiBc - I) REit)dN^i{t) 

a 

+ Y.^E{t) {M^2Bi - I) dN^2{t). (23) 



The Poisson increments dNai, satisfy Eqs. p4p and (jisp . 
and the transition rates are given by 



r 



1 



and 



F - Vr 

a 

r = Ti+Fa. 



(24) 

(25) 
(26) 



The quantities Lau and Mai, are real and positive func- 
tional of ipi, ■02 and Re- One has a great freedom in the 
choice of these functionals, the only restriction being the 
positivity. A definite choice will be made in the example 
below. 
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We observe that again the ip^, follow a pure jump pro- 
cess. If dNau{t) = 1 for a particular pair of indices a 
and V, which happens with probability Ta^dt, the state 
vector tpi, undergoes the jump 

~iLa^Aa^p„. (27) 



The Pauli spin operator of the central spin, constitut- 
ing the open system, is denoted by a with corresponding 
raising and lowering operators a± — ^((Ti ± ia2). The 
central spin couples to N environment spins with Pauli 
spin operators i — 1,2, . . . , N , through the raising 
and lowering operators 



At the same time Re carries out the jump 

Re MaiBaRs (28) 

if I' = 1 and the jump 

Re REM^2Bi (29) 

if = 2. Thus, for v = 1 the operator Ba acts from 
the left on Re, while for v — 2 the adjoint operator BI^ 
acts from the right on Re- Under the condition that all 
Poisson increments vanish, which occurs with probability 



1 - Taudt = 1 - Tdt, 



(30) 



we have d^p^{t) = and dRE{t) — TRE{t)dt, i.e., the 
"01, are left unchanged during dt while the environment 
matrix RE{t) follows a linear drift. 

Like in the case of the process constructed in Sec. Ill Bl it 
is easy to design an appropriate Monte Carlo algorithm 
for the stochastic differential equations and 
We note that in both cases the random matrix R{t) has 
the structure of a tensor product, which considerably re- 
duces the complexity of the problem. As a consequence 
of Eqs. pTjl and (|2ip the environmental states enter the 
expectation value for the reduced density matrix ps{t) 
only through the scalar product (X2(i)|xi(i)) or through 
the trace tr e Re (t). To simulate the non-Markovian dy- 
namics of an open system with these algorithms it thus 
suffices to record during the simulation of the process 
the various jumps and their moments of occurrence. At 
any time t the scalar product (X2(i)lxi(^)) the trace 
tTERsit) can then be expressed in terms of certain corre- 
lation functions of the environmental operators Ba- For 
many system-environment models the latter are known 
explicitly. This fact greatly facilitates the numerical im- 
plementation of the stochastic method. An example is 
discussed in Sec. IIII Cl (see, in particular, Eq. ([64])). 



III. APPLICATIONS 
A. The spin star model 

As a simple but instructive example of the Monte Carlo 
method, we investigate in the following a spin star model 
described by the time independent interaction Hamilto- 



H 



2A 



{<7+J^ + CJ-J+) . 



(31) 



N N 



(32) 



of the total angular momentum J of the environment. 
The initial state of the total system at time t = is 
taken to be a product state p{Q) = ps{0) Pe{0), where 
the reduced density matrix ps(0) of the central spin may 
be an arbitrary, possibly mixed state. The spin bath 
is assumed to be in an unpolarized infinite temperature 
initial state 



Pe{0) 



2-^/s, 



(33) 



where Ie denotes the unit matrix in He- 

This model can be solved analytically [13] . We express 
the solution for the reduced density matrix in terms of 
the components of the Bloch vector 

v{t)^tr{<7p{t)} (34) 

which are related to the reduced density matrix by 



Ps{t) 



P++it) P+-{t) 

l + V3{t) Vi{t)-iV2{t) 
Vl{t)+iV2{t) l~V3{t) 



(35) 



The components V3{t) and v±{t) = ^{vi ± iv2) are then 
given by the explicit expressions: 



v±iO) 
where 

P{j, m 

and 



= ^P(j,m) cos [2r(j,m)i], 



(36) 



P{h m) cos [V{j, m)t\ cos [T{j, ~m)t] , (37) 



1 

2^ 



N 



N 



+ J 



N 



N 

+J + 1 



r(j, m) = 2A 



iii + 1) ^ m(m + 1) 



N 



(38) 



(39) 



These expressions may be obtained by solving the 
Schrodinger equation of the model with the help of the 
fact that the manifolds spanned by the states |-|-) (g) |j, to) 
and |— |j, TO-l-1) are invariant under the time evolution. 
Here, |±) are the eigenstates of 173 with corresponding 
eigenvalues ±1, and \j,m) denotes an eigenstate of the 
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square J ^ of the angular momentum J of the spin bath 
and of its 3-component J3 with respective eigenvalues 
j{j + 1) and m. 

The quantity P{j, m) defined in Eq. ([38|) is the proba- 
bility of finding the quantum numbers j and m in the ini- 
tial mixture representing the state ((33l) p^ . As usual, for 
N even j takes on the values j = 0, 1, 2, . . . , N/2, and the 
values j = f/2,3/2,...,iV/2ifA^isodd. For a given j the 
quantum number m takes the values — j, — j -f 1, . . . , + j. 
It is easy to check that the probability distribution 
P{j, m) is normalized as follows, 



^P(j,m) = 1. 



(40) 



In the limit — > 00 of an infinite number of bath spins 
the above formulae lead to the asymptotic expressions 



l + 2g(i). 



(41) 
(42) 



where 



with X = \plAt. The function erfi(a::) denotes the imag- 
inary error function, which is a real-valued function de- 
fined by 



erfi(a;) 



_ erf(zx) _ J_ 



I V^^A:!(2fc + l)' 



(44) 



B. The representation i? = I $1) ($2 1 

We first illustrate the stochastic representation for the 
process defined in Sec. Ill Bl To this end, the initial state 
psp is realized through a mixture of the states |j, to), 
where the quantum numbers j and m follow the joint 
probability distribution P{j,m) given in Eq. (|38p . To 
analyze the process we therefore have to describe the 
stochastic evolution of the initial states 



1*1(0)) = |$2(0)> = |+>®|i,m). 



1*1(0)) = |$2(0)) = |-)®b-,m), 



(45) 



(46) 



from which we can reconstruct the reduced density ma- 
trix of the central spin. 

Let us consider first the initial state (US)) . According to 
the interaction Hamiltonian (PT]) the index a in Eq. (|3|) 



assumes two values a = ± with corresponding time in- 
dependent operators 



B± = 



2A 



J: 



(47) 



The jumps of the process thus take the form 

llx." 



lk±V'.|| 



•hXu- (48) 



Since ipv{0) = |+) the states ^„ jump between the states 
1+) and |— ), whereby each jump contributes an addi- 
tional factor of (— j). Let us denote the number of jumps 
of |<i>,^) during the time interval from to f by = n^{t). 
We then have ipu{t) ~ (~*)""l+) if "-i^ is even, and 
ipvit) = if f^v is odd. The corresponding envi- 

ronment states are given by 



for even Ui,, and by 

X.(i)-|j, m + l)ei^(^''")* 



(49) 



(50) 



for odd n,j. The rates T± of the jumps are determined 
as follows. 



(51) 



2A 



r+ = -=||J_|j,m + l)||=r(j,m). (52) 



Thus we have r_|_ = r_ = r(j, to), where r(j, m) is given 
by Eq. ([39]) . The deterministic drift of the process there- 
fore yields a factor exp[r(j, m)t\, which has already been 
taken into account in Eqs. (|49)) and (jSO]) . 

On using this information we can determine the dy- 
namics of the populations of the reduced density ma- 
trix. Considering psiO) = we have by virtue of 
Eq. dni): 

p++{t) = {+\Ps{t)\+) (53) 

= Tii{+\Mt))mt)\+){x2{t)\xim 

= E(^^;(nl,r^2)(-^)"l(+^)"'e2^0^™)*) . 

Of course, a given realization of the process contributes 
to the expectation value only if ni and 712 are even. This 
fact is accounted for by the first factor w(ni, 712) which is 
defined to be equal to 1 if ni and 712 are even, and equal 
to zero otherwise. The second and the third factor un- 
der the expectation value take into account the jumps of 
(factor (-i)"0 and of V2 (factor [(-i)"']* = (+«)"')• 
Finally, the exponential function represents the contri- 
butions from the scalar product (X2(i)lxi(0) which, ac- 
cording to Eq. equals exp[2r(j, to)<]. 

It is clear from the general theory outlined in Sec. Ill Bl 
that Eq. (|53p is an exact representation of the popula- 
tions. Nevertheless, it might be instructive to see ex- 
plicitly how the exact solution p6p for the 3-component 
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of the Bloch vector emerges from the expectation value 
([55)1 . To this end, we note that the states evolve 
independently and that the transition rates of the pro- 
cess are time independent. This implies that the random 
numbers ni(t) and 712 (t) follow independent Poisson dis- 
tributions with the same mean value of T{j, m)t: 



(54) 



The expectation value therefore becomes 



xe2r(^'''")*P(ni,i)P(n2,t). (55) 

The first sum extends over all possible values of the 
quantum numbers j and m occurring in the initial state 
(see Sec. IIII Al) . while the second sum runs over all 
'ni,n2 = 0, 2,4,.. .. Substituting the expression ([54)1 into 
Eq. ([55]) . we get 

p++it) = ^F(j, m) ^(-z)"n+0"^ 



[r{j,m)tr- [r{j,m)tY 



ni'. 



712! 



^P(j,m) cos2(r(j,m)i). 



(56) 



Using, finally, the relation V3{t) — 2/9++ — 1 we see that 
Eq. ([56)1 leads to the exact expression p6|) for the 3- 
component of the Bloch vector. Thus we see explicitly 
that the stochastic process indeed reproduces correctly 
the exact time evolution of the system. 

In a similar way one finds the coherence f _ (t) = (t) 

of the central spin. To this end, we have to consider 
also the initial state . The resulting process is essen- 
tially the same as above, with the only difference that 
now r+ = r_ = r(j, — to). To find the expectation 
value representing v_ (t) we have to use the initial states 
|$i(0)) = \+) (g) \j,m) and |$2(0)) = |-) (S \j,m). With 
the initial condition v- (0) = 1 we then have 



v^{t) = {+\ps{t)\-) 

= E{{+\Mt))mt)\-){x2{t)\xim 



(57) 



= E w 



7(ni, n2)(-i)"' (+i)"^e(^(j'^™)+^(j'^-"'»*) 



It is easy to verify that this expectation value leads to 
the exact expression (|57|) for the coherence of the central 
spin. 

Figure[T]shows an example of a Monte Carlo simulation 
of the stochastic process defined by the differential equa- 
tions and (|14p . As can be seen, the Monte Carlo 
simulation reproduces the exact solution with high ac- 
curacy over the range of time shown. The figure also 
indicates the growth of the size of the statistical errors 
which have been estimated from the sample of realiza- 
tions generated. 



0.5 



-0.5 



0.2 



0.4 



At 



0.6 



0.8 



FIG. 1: Monte Carlo simulation of Eqs. (fT2|l and (IT4)l using a 
sample of 10^ realizations: The 3-component V3 = 2p++ — 1 
(dots and errorbars) of the Bloch vector and the coherence 
V- = p+- (crosses and errorbars) of the central spin coupled 
to a bath of 10^ spins through the Hamiltonian (|31|l . The con- 
tinuous line and the broken line show the analytical solution 
given by Eqs. (|36|) and p7|l . respectively. 



Beyond the point At k, 1 the statistical errors strongly 
increase. This is a typical feature of the Monte Carlo sim- 
ulation method which also appears in many other models. 
In fact, if one measures the size of the statistical errors by 
means of the Hilbert-Schmidt distance between the ran- 
dom operator R{t) and its mean value p{t), one can show 
9] that for large times the fluctuations grow roughly as 
exp(2roi), where To represents an upper bound for the 
rates T^. As can be seen from the above example this 
exponential increase of the errors is mainly due to the 
corresponding increase of the norm of the environmental 
states Xi^- 



C. The representation R = \'4>\) {'ii)2\ ® Re 



Let us now analyze the process defined in Sec. Ill C[ 
which is particularly suited to simulate the limit of an 
infinite number N of bath spins. We choose the quanti- 



\\Ac.M\' 



(58) 



In our example we then have Lau = 1- Thus, V'l^ again 
jumps between states |-|-) and |— ), whereby each jump 
contributes a factor of (— i). Hence, we have 

p++(t) = E(^«(nl,r^2)(-^)"^(-l-^)"^tr^^i?B(^)). (59) 

The aim is now to determine the trace over the random 
environmental operator Rsit). This will be done in the 
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limit N oo. In this limit we have for a fixed number 



{{j+j-f) = k\ 



where we define 

(0> =trB{Op£;(0)} = 2-^tr£;O 



(60) 



(61) 



for any bath operator O. We further choose the jump 
rates Fi = r2 = y/2A (the factor \f2 is introduced for 
convenience), such that F = Fi + r2 = 2^/2A. The drift 
contribution to Rsit) is therefore given by exp{2\/2 At). 
The jumps of the random matrix Re take the form 



V2A 



B±Re 



JzrRf 



^ReB^-^-ReJ,. 



(62) 



(63) 



Since Re{0) is proportional to the identity and since the 
order of application of the operators J± is irrelevant in 
the limit N oo, we conclude that 



trERsit) 




^] {{J+J-f)e^\ (64) 



where we have defined k ~ (jii + n2)/2, assuming that 
both ni and n2 are even. Employing Eq. (|60p we there- 
fore get 

tTERsit) = fc!e". (65) 

Hence, the expectation value ((59|) becomes: 

p++{t) = E (?«(ni,n2)(-l)^fc!e") . (66) 

It is again instructive to see explicitly how this expres- 
sion leads to the formulae (|¥T|) and Using the fact 
that ni(t) and n2{t) are independent and follow Poisson 
distributions with mean value ^/2At, one finds 



E 

rii ,712 

oo 

E 

k=0 



(V2At)"i (V2Ai)"2 



712! 



(-l)^-fc! 



"E 



We recall that the sum in the first line extends over the 
values ni, 712 = 0, 2, 4, . . ., and that we use the definition 
k = (ill +7i2)/2. The second sum in the second line runs 
over the values ni = 0, 2, . . . , 2k. This sum is found to be 
equal to 1 for fc = and equal to 2^'^"^ for fc = 1, 2, 3, . . .. 
Thus we obtain 



p++{t) = 1 + 



E 

k=l 



(-l)fefc! 
2(2fc)! 



i2V2At)'^'' = l + git), (67) 




FIG. 2: The 3-component 113 = 2p++ — 1 of the Bloch vector of 
the central spin coupled to an infinite number of bath spins 
through the Hamiltonian (|31|) . Dots and errorbars: Monte 
Carlo simulation based on the stochastic differential equa- 
tions (|22|) and (|23p with 10^ realizations. Continuous line: 
Analytical solution given by Eqs. (|41|l and (|43|l . 



from which we get V3{t) = 2p++{t) — 1 = 1 + 2g{t), 
where the function g{t) has been introduced in Eq. (|43p . 
A similar reasoning leads to the relation v±{t) = 1 -f 
g{t). These results coincide with those obtained from 
the solution of the Schrodinger equation [see Sec. IIII A] . 

Figure O shows the results of a Monte Carlo simula- 
tion of the stochastic process defined by the differential 
equations (|22p and The Monte Carlo simulation re- 
produces the exact solution with high accuracy over the 
range of time shown, and we again observe the growth of 
the statistical errors. 



IV. CONCLUSIONS 

We have investigated two methods that yield an ex- 
act stochastic unravelling of the non-Markovian quantum 
dynamics of open systems by means of a piecewise deter- 
ministic Markov process. These methods yield Monte 
Carlo simulation techniques that are generally applica- 
ble for the investigation of the short-time behavior of 
the open system's dynamics. Due to a possible expo- 
nential increase of statistical fluctuations for large times 
a numerical simulation of the long time behavior is, in 
general, impossible in practice. 

However, a great advantage of the method is given by 
the fact that it is exact and that it allows the treatment of 
arbitrary correlations in the initial state. It must be em- 
phasized that a Monte Carlo simulation not only yields 
an estimate for the desired averages, but also for the sta- 
tistical errors. As long as the latter are small the tech- 
nique yields excellent predictions about the short-time 
behavior and thus offers the possibility to control and as- 
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sess the performance of other approaches and approxima- 
tion schemes. In particular, the method may find impor- 
tant apphcations in the simulation of non-Markovian de- 
coherence phenomena which are dominated by the short 
time behavior of the open system. 

The formulation of the stochastic simulation method 
has been given here in the interaction picture, assuming 
that the free dynamics of the system and the environ- 
ment are known. If this is not the case one can use an 
analogous formulation of the stochastic dynamics in the 
Schrodinger picture which includes the given Hamilto- 
nian operators for the system and the environment into 
the deterministic drift of the stochastic differential equa- 
tions [§]. 

In our examples we have restricted ourselves to the 
case of an unpolarized (infinite temperature) initial state 
of the spin bath. It should be noted that the stochastic 
method is also applicable to polarized initial states. For 
instance, one can consider an initial equilibrium state of 
finite temperature. Introducing a spin bath Hamiltonian 
of the form He — loJ^, one then has to multiply the 
probability distribution P{j, m) defined in Eq. (|38|) with 
the m-dependent Boltzmann factor exp(— /Jwm), where 
f3 is the inverse temperature. 

There are two basic strategies for the improvement of 
the Monte Carlo technique. The first one employs the 
freedom in the choice of the stochastic time-evolution in 
order to minimize the statistical errors [ll|, This 



can be done by an appropriate modification of the noise 
terms of the stochastic differential equations. A further 
possibility is to introduce additional terms in the deter- 
ministic part of the equations of motion. This approach 
leads to a stochastic mean field dynamics which is sim- 
ilar to the one used in the Monte Carlo wave function 
method for interacting many-body systems [l5l | . 

The second strategy is to reduce the size of the sta- 
tistical fluctuations by using more complicated stochas- 
tic functionals whose expectation values lead to the re- 
duced density matrix. The methods investigated here 
represent the correlated states of the composite system 
through the average over random operators with a spe- 
cific given structure, namely a tensor product structure. 
This ansatz does by no means exhaust all possibilities. 
There are many other possible ways of constructing an 
exact stochastic representation of the dynamics which 
seem worth being explored in a more systematic manner. 
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